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ABSTRACT 

This paper constructs an analytic form for a triaxial potential that describes the 
dynamics of a wide variety of astrophysical systems, including the inner portions of dark 
matter halos, the central regions of galactic bulges, and young embedded star clusters. 
Specifically, this potential results from a density profile of the form p{m) oc m"^, where 
the radial coordinate is generalized to triaxial form so that m? = ja^ + /b'^ + /(? . 
Using the resulting analytic form of the potential, and the corresponding force laws, 
we construct orbit solutions and show that a robust orbit instability exists in these 
systems. For orbits initially confined to any of the three principal planes, the motion 
in the perpendicular direction can be unstable. We discuss the range of parameter 
space for which these orbits are unstable, find the growth rates and saturation levels 
of the instability, and develop a set of analytic model equations that elucidate the 
essential physics of the instability mechanism. This orbit instability has a large number 
of astrophysical implications and applications, including understanding the formation 
of dark matter halos, the structure of galactic bulges, the survival of tidal streams, and 
the early evolution of embedded star clusters. 

Subject headings: stellar dynamics, celestial mechanics — methods: analytical — galax- 
ies: kinematics and dynamics, halos — stars: formation 



1. INTRODUCTION 

Many types of astrophysical objects are essentially collisionless systems with extended mass 
distributions, including dark matter halos, elliptical galaxies, galactic bulges, and star clusters. Al- 
though many of these systems display triaxial shapes, explicit analytic models for the gravitational 
potentials, force laws, and orbits of the triaxial incarnations of these systems are rare. This paper 
constructs an analytic form for the triaxial potential of an important class of systems, namely those 
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in which the density profile approaches the form p '-^ 1/r in the spherical limit. With the potential 
in hand, we find the corresponding force terms (see also Poon & Merritt 2001; hereafter PMOl) 
and study the orbits allowed by this triaxial system. One contribution of this present work is the 
discovery of a robust instability, wherein the orbital motion initially confined to any of the three 
principal planes can be unstable to the growth of the amplitude of the motion in the perpendicular 
direction. In addition to establishing the existence of the instability, this paper develops a set of 
model equations that allow for a (mostly) analytic description of the instability, the magnitude of 
the growth rate, and the saturation mechanism. 

This particular form for the density profile p ~ 1/r arises in many contexts. The Hernquist 
density profile, which has this form in the inner limit (Hernquist 1990), was originally constructed 
to describe the i?^/^ law in galactic bulges and elliptical galaxies. Since the Hernquist profile arises 
in many other contexts, we previously studied orbital solutions in the full (spherical) Hernquist 
potential and obtained a number of analytic results that are applicable to the present analysis 
(Adams & Bloch 2005, hereafter AB05). However, these bulge systems are known to be triaxial 
(Binney k, Trcmaine 1987, hereafter BT87; Binncy & Merrifield 1998), and hence the triaxial 
generalization of the potential is important for understanding these galactic-scale systems. For 
example, triaxial potentials support box orbits, which allow stars to wander arbitrarily close to the 
galactic center; this behavior affects the rate at which stars interact with the central black hole. 

On a larger scale, numerical simulations indicate that the density profiles of dark matter halos 
display a nearly universal form in which p ~ 1/r in the inner regime (Navarro et al. 1997, hereafter 
NFW). Recent numerical simulations that carry the calculations into the future (Busha et al. 2005) 
show that the asymptotic form of the density profile becomes steeper at large radii and is well- 
described by a Hernquist profile (Hernquist 1990) but maintains the same p ~ 1/r form on the 
inside. Additional studies suggest that this result for the inner form is robust; if the density profile is 
written in the form p oc r^^, then p < 1.2 (see Diemand ct al. 2004, Hayashi & Navarro 2006, Huss 
et al. 1999, and references therein). Nonetheless, these systems are triaxial, rather than perfectly 
spherically symmetric, with typical axis ratios of 5:4:3 (e.g., Kasun Sz Evrard 2005, Jing & Suto 
2002). A host of other authors have considered halo shapes, and always find roughly ellipsoidal 
shapes, but considerable variation exists; in addition, the axis ratios can vary with radius within the 
halo (sec Allgood et al. 2006 for further discussion). Another result from numerical simulations is 
that the inner portions of dark matter halos tend to have (mostly) isotropic velocity distributions, 
while the outer regions display more radial distributions; the orbit instability considered here acts 
to isotropize the inner regions of these systems and may help explain this finding. 

On a smaller scale, the density profiles of cluster forming molecular cloud cores are observed 
to have substantial non-thermal line-widths (Larson 1985, Jijina et al. 1999); the size versus line- 
width relation can be used to infer an effective equation of state, which, in turn, implies a density 
profile of the form p ^ 1/r in the spherical limit (see Lizano & Shu 1989; Jijina & Adams 1996; 
McLaughlin & Pudritz 1996). These cores often form small star clusters (N = 30 — 3000) and the 
stellar members orbit within the corresponding potential, which is dominated by the remaining gas 
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(that has not turned into stars) . These orbits determine the interaction rate between young stellar 

objects and also determine their radiation exposure (both processes generally have a destructive 
effect - see, e.g., Adams ct al. 2006). These systems are also observed to be triaxial, and roughly 
spheroidal, with typical axis ratios of order 2:1 (Myers et al. 1991, Ryden 1996). In addition, recent 
observations suggest that some systems are highly flattened, with even more extreme axis ratios 
(P. Myers, private communication). The triaxial nature of the potential allows star/disk systems to 
execute box orbits, which brings them close to the cluster center where massive stars produce large 
amounts of radiation (which can destroy the disks). On the other hand, orbit instability acts to 
round out the clusters, even in the absence of stellar scattering encounters. In any case, a complete 
description of the dynamics of young embedded star clusters requires an understanding of orbits in 
these triaxial systems. 

Orbits and orbital instabilities in extended collisionless systems have a long history. As is 
well known (BT87), triaxial systems display a wide variety of orbits and a number of authors 
have studied orbits in systems similar to those considered here (e.g., de Zeeuw & Merritt 1983, 
Statler 1987; Holley-Bockelmann et al. 2001, 2002; PMOl, Poon k Merritt 2002, Terzic & Sprague 
2007, and many others). The radial orbit instability, where nearly radial orbits are unstable to 
perpendicular motions, was first observed in numerical experiments by Henon (1973). This orbit 
instability was considered further by several authors (including Barnes et al. 1986, Aguilar &: 
Merritt 1985, MacMillan et al. 2006). The instability of motion in the directions perpendicular 
to the principal planes was considered by Binney (1981); subsequent work shows that this type 
of behavior arises in triaxial potentials such as those considered herein (e.g., Merritt &; Pridman 
1996). However, a detailed (and analytic) description of the instability mechanism has not been 
given and this issue represents a main focus of this paper. 

Study of instability in extended mass distributions can be carried out in several ways: In one 
case, the global stability of the system is considered and perturbation theory is used to study the 
evolution of the distribution function (e.g., Aguilar & Merritt 1985; Barnes et al. 1986; Palmer &; 
Papaloizou 1987). Instead of considering the distribution function, another approach is to consider 
the potential as fixed and study the perturbations of the orbits themselves (e.g., Binney 1981, 
Scufiaire 1995, Papaphilippou & Laskar 1998). Note that these two types of instability are related, 
but are not equivalent. If the distribution function is unstable, as in the former case, then the 
orbits must change (and hence be unstable in a sense); it is possible, however, for particular orbits 
to be unstable while the distribution function remains stable. In the context of individual orbits, 
an important issue is to study the degree of chaos of particular orbits, generally by computing 
the lyapunov exponents (e.g., PMOl, Valluri &; Merritt 1998). Although one expects more chaotic 
orbits to lead to greater "instability", the two issues are not equivalent. For example, orbits that 
are confined to a given plane can be chaotic and even have large lyapunov exponents, but could 
nonetheless remain stable to perturbations in the perpendicular direction. The study of instability 
of individual orbits thus represents a separate line of inquiry. Here we adopt the latter approach, 
follow the instability of individual orbits, and develop a detailed description of the instability 
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mechanism. 

This paper builds on the aforementioned previous studies (see also de Zeeuw & Pfenniger 
1988), and presents new results concerning the instability of orbits in triaxial potentials. We note 
that a great deal of the previous work is either numerical or uses a simple logarithmic model 
for the potential, where <^ oc ln[a;^/a^ + + z^/c^]. Although the numerical work cleanly 

demonstrates the nature of the orbits and the existence of instability, it does not elucidate the 
detailed mechanisms that lead to the instability. Previous analytic work using logarithmic potential 
is limited to the study of axis ratios close to unity because the corresponding density (e.g., obtained 
from = AttGp) becomes negative for more extreme axis ratios (although more complicated 
generalizations allow for more extreme axis ratios - see Schwarzschild 1993). In particular, the 
estimated values of a : 6 : c = 5:4:3 for dark matter halos (Kasun & Evrard 2005) are over the limit 
for the applicability of the usual logarithmic model potential. For example, if we scale the axes so 
that a = 1, the requirement for the logarithmic potential to have nonnegative density everywhere 
can be written in the form c > 6/vT+~P^. For b = 4/5, typical of dark matter halos, this constraint 
becomes c > 4/\/4T 0.62 (i.e., the typical the axis ratios found in numerical simulations just 
fail to satisfy this constraint). For b = 1/2, typical of cluster forming molecular cloud cores, the 
constraint is more restrictive, c > l/\/5 ~ 0.45 (just below b = 1/2), whereas we expect even more 
extreme axis ratios (perhaps as small as c ~ 0.10). 

Although the results of this work are applicable to a wide range of astronomical problems 
(outlined above), this paper focuses on the construction of the triaxial potential, a brief description 
of its orbits, and a detailed analysis of the orbit instability. The paper is organized as follows. 
In §2 we construct the analytic form for the potential resulting from the triaxial generalization of 
a density profile p ~ 1/r, and present the corresponding analytic expressions for the forces and 
axis ratios. In §3 we briefly survey the possible orbit solutions in the triaxial potential. Next we 
study the stability of orbits, and show that orbits in all three principal planes can be unstable to 
motion in the perpendicular direction (§4); we also develop analytic model equations to describe the 
instability mechanism (§5). The paper concludes, in §6, with a summary of results and a discussion 
of their astrophysical implications. Although this paper focuses on the inner limit (where p ~ ^/r), 
Appendix A presents analytic expressions for the potential and force terms valid over the full radial 
extent of the axisymmetric version of the NFW profile. 



2. ANALYTIC POTENTIAL AND FORCES 



The overarching goal of this study is to understand orbits in the potentials resulting from a 
generalized density profile of the form 

fim) , . 

Pgen = PO , (1) 
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where po is a density scale and where the variable m has a triaxial form 

2 x2 y2 2 

0^ 

Keep in mind that the variable m defines the surfaces of constant density, whereas the radial 
coordinate ^ is given by 

f = x'^ + y'^ + z\ (3) 

The function f{m) is assumed to approach unity as m ^ so that the density profile has the 
form p ~ 1/m in this inner limit. Two standard forms for extended mass distributions (with many 
astrophysical applications) are the NFW profile (Navarro et al. 1997) and the Hernquist profile 
(Hernquist 1990), which have density profiles of the forms 

Pnfw = —7— and phq = —7— ^3 • (4) 

For the rest of this analysis, we use a dimensionless formulation so that po = 1. In the spherical 
limit, m = r/r^, where is a scale radius (NFW). In this treatment, we normalize all length scales 
by rg to make the variables {x,y,z) and the geometric constants {a,b,c) dimensionless. 

For any density profile of this general (triaxial) form, the potential can be found by evaluating 
the integral 

9{x,y,z)-Ij^ (a2 + n)i/2(c2 + ^)i/2(52 + ^)i/2 ' 

where m is given by equation ^ where the function 'ilj{m) is defined through the relation (BT87, 
Chandrasekhar 1969) 

POO 

ilj{m) = / 2mdmp{m) . (6) 



We then define a constant A by 

A=2 I f{m)dm, (7) 







where we assume that the integral converges, so that the function ip takes the form ip = A — 2m in 
the inner limit (m — > 0). For the NFW (Hernquist) profile, the integral is easily evaluated and the 
constant ^ = 2 (1). Keep in mind that in this dimensionless formulation, we have set a number of 
constants to unity and let the potential have a positive value. 



2.1. Axisymmetric Inner Limiting Form 

For the axisymmetric problem, we can set a = 1 = b and define = + After defining 
two new "constants" according to 

7 = 1 - and q = jR^/f , (8) 
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the potential in the inner limit can be written in the form 



cos 



2z, 
H In 

7 



4e, 



7 



1), 



(9) 



where the constant A depends on the form of the density profile over the rest of its range (as defined 
above). We note that for the axisymmetric version of the NFW profile, the full potential (valid for 
all m, not just in the inner limit) can be found analytically. This result is presented in Appendix 
A, and equation ([9]) agrees with the inner limit of the general NFW potential. 



2.2. Triaxial Inner Limiting Form 

Next we generalize to the triaxial case, where a > b > c > 0. If we consider the inner limit of 
a density profile of the form of equation ([1]), then the potential can be written as two terms 

^ = Ah- 2I2 , (10) 

where A has the same meaning as before and where 

h ^ r , (11) 

and 

[eV + An + r]^/2 



Jo = / du 7-77; ——^ 7 . (12) 

Jq (a2+u)(62+u)(c2+n) ^ ^ 

The radial coordinate is defined through equation ([3]) and the remaining coefficients in the nu- 
merator are given by the following functions of the coordinates, 

A = (&2 + c2)x2 + (a^ + c2)y2 + (a2 + ^2)^2 ^ ^ ^2 ^2^2 ^ ^2^2y2 ^ ^2^2^2 ^ (-^3) 

Note that the first integral Ii defines the depth of the potential well. Further, Ii does not 
depend on the spatial coordinates — it is a constant for a given set of axis ratios (a, b, c) — and is 
determined by the obliquity parameter rjoh = (a^ — 6^)/(a^ — c^). 

The second integral I2 can be expanded into three terms according to 

1 /•°° , [c2^2^Au + r]^/2 

' au- 



(a2 — 52) (^2 _ ^2) a2 + n 



' au- 



(a2 -62)(62 _c2) 7o b'^ + u 

+ (a2 _ c2)(52 _ c2) c2+n • ^ ' 
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The integrals can be evaluated individually and recombined to take the form 



F{b)ia' 



(a2-62)(a2_c2)(52_ 
Ab^ - 2^%^ + 2F{b)Cb^ 




2r - Aa2 + 2F{a)Vf 
Ac2 - 2^2c4 + 2F(c)ec2 



+ 



2r - A62 + 2F(6)Vr 
where the function F is defined by 

F{a) = [^a^-Aa^ + rf^ 



2T - Ac2 + 2F{c)Vf 



(15) 



(16) 



The spatial dependence of the potential is contained in the functions ^, A, and T, which depend on 
the usual coordinates (x, y, z) through the relations of equation (fT3]) . Notice that the second term, 
as written, has both real and imaginary parts, since F'^{b) is negative (under the usual ordering 
a> b > c). We get agreement between the analytic form and the numerically evaluated form when 
we take the real part of equation (jlSp . Alternatively, we can rewrite the second term so that the 
potential integral takes the form 



(a2-62)(a2_c2)(52_c5 

+\F{b)\{a^ - c^) sin-i 



F{a){c^ 



b^)\n 



2c2 



A - 2¥i 



VA2 - 4e2r 

+ F{c){b^ - a^)\n 



Aa2 - 2g2Q4 + 2F(a)ga2 
2r-Aa2 + 2F(a)Vf 

. 2r/fe2-A ^ 

- sm = 

Ac2 - 2i^c^ + 2F{c)i^ 



2T - Ac2 + 2F{c)^/f 



(17) 



where the second term is now manifestly real. We also note that in the evaluation of the second 
term, it is sometimes advantageous to use a trigonometric identity to write the sin~^ terms as tan~^ 
expressions. 



2.3. Force Terms 

The components of the force can be obtained by direct differentiation of the potential (eq. 
|17|). Alternatively, one can start with the original integral expression for the potential, differentiate 
first, and then perform the integration. The second procedure is simpler, but both result in force 
components of the following forms (where we have re- introduced the minus sign): 



) £ 

dx 



2x 



■In 



2F(a)^/r + 2r- Aa^ 



[2F{a)i + A - 202^2 



(18) 



,dl2 
' dy 



2y 



\m\ 



sm I 



,2c2 



26^e 



VA2 - 4^2? 



sm 



-1 



2r/62-A X 

VA2 -4e2r/ 



(19) 
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dh 

az 



2z 
F(c) 



In 



2F{c)VT + 2T - Ac"^ 



c2[2F(c)e + A-2c2e2 



(20) 



We note that equivalent expressions for the force terms have been derived previously (PMOl), and 
that the results can be expressed in a variety of forms. In particular, the sin~^ functions can 
be written as tan^^ functions; for example, in ther first term of equation (|19p. @ = sin~^[(A — 



2^^C^)/\/A2 - 4^2r] = tan-^[(A - 252^2)/2^VA62 - T - 62]. Both forms are useful for numerical 
evaluation of the forces, depending on the context. In addition, the leading coefficients of the force 
terms obey identities that allow for different expressions, e.g., 2x/F{a) = 2{a^ — b'^)~'^/'^{a?' — (?)~^/'^ ^ 
with similar forms for the other forces. 



2.4. Shape of Equipotential Contours 



One common feature of triaxial systems is that the surfaces of constant potential are generally 
rounder (closer to spherical symmetry) than the surfaces of constant density, which are ellipsoidal 
(by construction) for the class of systems considered here. In these systems, the surfaces of constant 
potential are given by setting equation (jl7p equal to a constant; although this form is analytic, it 
remains both implicit and algebraically complicated. However, we can find the "axis ratios" for the 
potential surfaces by evaluating the potential along each of the principal axes. For example, along 
the X-axis, the potential reduces to the form 



$ = $o-xJ^(a,6,c), 
where all of the geometry is encapsulated in the function 



J cr. 



(a2 _ c2)l/2(a2 - 62)1/2 



In 



6(a2 



r.2U/2 



c a 



62)1/2 



a[{a? - c^fl-^ - (o2 - 62)1/2] 



(21) 



(22) 



Along the other two axes, the potential takes a similar form with 



and 



62)1/2(62 



U/2 { 



sm 



a2 +c2 



262 



sm 



a2(c2_52)^^2(^ 



}, (23) 



J. 



-2 



(a2 _c2)l/2(52 _c2)l/2 



In 



6(^2 _ ^2)1/2 _ ^(52 _ ^2)1/2 



c[(a2 



c2)l/2 _ (52 _ ^2)1/2] 



(24) 



The surfaces of constant density are ellipsoids. For example, if we want to find the locations 
on the X-axis and the z-axis with the same density, we require xja = zjc or zjx = c/a (where 
we take positive values). Proceeding in analogous fashion, if we want to find the locations on 
the X-axis and the z-axis with the same value of the potential, we require xJr^ = zJz so that 
z/x = Jx/Jz- The other axis ratios can be found similarly. Further, the above equations imply 
that c/a < Jx{a,b,c)/ Jz{a,b,c), and similarly for the other pairs of axes, so that the surfaces of 



- 9 - 



constant potential are indeed rounder (closer to spherical) than the surfaces of constant density. 
Moreover, we have obtained analytic expressions for the axis ratios. 

In the limit of extreme axis ratios, the above expressions simplify considerably. In the limit 
a S> 1 and c ^ 1, the integrals that define the shape of the potential approach the asymptotic 
forms 

J^^'-MM, J,^l, and ._,^1M?Z£). (25) 

a a 

In the opposite limit, where a — > 1 and c — > 1, the potential becomes spherical. 



2.5. Nearly Spherical Density Profiles 

In some applications, it is useful to have the form of the potential for a density distribution 
that is nearly spherical, but nonetheless displays triaxial departures from perfect symmetry. For 
the sake of definiteness, this subsection considers the triaxial density profile with axis ratios 

= l + 5a b = l, and = 1 - Sc , (26) 

where physically meaningful solutions require < 6a, 6c < 1, but this form is most useful in the limit 
6a,Sc ^ ^- After propagating this ansatz through the equations developed above, the potential can 
be written in the relatively simple form 

u 2 8^ ' ^ ^ 

where ^ is the (usual) radial coordinate given by equation ([3]). The density corresponding to this 
potential takes the form 

1 + 2 {6c - 6a) {6aX^ - 6cZ^ 

It is straightforward to verify that this form is the leading order correction to the general density 
profile of equation ([1]) in the limit of small departures from sphericity. Further, this density field will 
be positive as long as 6a < 1/2 + 35c/4, i.e., for sufficiently spherical systems. For more flattened 
profiles, the full form of the equations presented in the previous sections must be used. Note that 
this potential-density pair has a density profile of the form p ~ 1/^ in the spherical limit, rather 
than the more commonly used logarithmic potential, which corresponds to a density profile of the 
form /3 ~ 1/^^. As a result, this form provides a useful alternative to the logarithmic potential and 
can be applied to many astrophysical systems of interest. Notice also that the density/potential 
pair found above simplifies even further for the symmetric case where 6c = 6a- 



P - 7- + • v2»j 



3. ORBITS 



A host of previous authors have studied the various orbits that are supported by triaxial 
potentials such as those considered herein (e.g., see BT87, Richstone 1982, Hunter & de Zeeuw 
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1992; Schwarzschild 1993, PMOl, Holley-Bockelmann et al. 2002). A wide variety of such orbits 
exist, including radial orbits, tube orbits, box orbits, and resonant orbits. A full accounting and 
presentation of all of the possible orbits has (mostly) been covered by the aforementioned previous 
work and is not the goal of this paper. Instead, we are primarily concerned with studying the orbit 
instability that arises in these systems, as well as obtaining analytic results whenever possible. 
This section develops analytic solutions for principal axis orbits (§3.1) and for radial orbits in the 
spherical limit (§3.2). A brief discussion of the surfaces of section, the fraction of box orbits versus 
loop orbits, and the dependence of these results on the axis ratios is presented in §3.3. 

Throughout this paper, when numerical integration of the orbits is required, we use a Bulirsch- 
Stoer (B-S) integration scheme (Press et al. 1992). This method is both accurate and explicit. For 
the systems at hand, our B-S scheme incurs errors in relative accuracy of order 1 part in 10^^ per 
total time step, with typical accumulated errors of 1 part in 10^, small enough that numerical error 
is not an issue. For completeness, we note that some longer runs accumulate a total error of 1 part 
in 10^. 



3.1. Principal Axis Orbits 



For orbits that are constrained to one of the principal axes, the equation of motion reduces to 
the simple form 



(fwj 



-Qj = constant , 



(29) 



where Wj represents the x, y, or z coordinate. The equation is thus the same as that of a baseball 
thrown vertically in the Earth's gravitational field (with no air resistance), although the trajectories 
repeat (oscillate) in the present application. The constant force strengths follow directly from 
equations (US -[20]). 



For orbits along the principal axes, the force terms are constant and can most easily be written 
in integral form 

(30) 



du 

(^/ + a2)[(n + 62)(n + c2)]^/^ ' 

/•°° du 

^^~Jo (n + 62)[(u + a2)(n + c2)]^/^ ' 

/■°° du 

^'~Jo ('u + c2)[(u + a2)(n + 52)]V2 • 

These integrals can be evaluated. For example, the force along the x-axis can be written 

1 



(31) 
(32) 



■In 



2[7x + bc^]/a'^ + 62 + c2 - 2a2 



2^ + 62 + 



2a2 



(33) 
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where we have defined = (o^ — 6^)(a^ — c^). As before, the integral along the intermediate axis 
takes a different form, 



sm 



a2 + ^2 _ 252 



cP- — (? 



sm 



+ - 26^ + 
— 



(34) 



where we define = (a — b ){b — c ), whereas the integral along the shortest axis becomes 



In 



2hz + abJYAIc +a' + b'- 2c 



2 VtI + a2 + 62 _ 2c2 



(35) 



where 7^ = (a^ - c^){¥ - c^). 



3.2. Radial Orbits in the Spherical Limit 

In the limit of a spherical Hernquist potential, the radial equation of motion takes the simple 

form 

dt2 +(1 + ^)2 W 

where ^ = r/r^ is the dimensionless radial coordinate and where 

Lol = 2ttGpo . (37) 

The parameter ojq thus sets the fundamental time unit for this potential. This equation of motion 
can be directly integrated to find the solution. For the case of an orbit starting with initial radial 
coordinate ^0 and zero (radial) velocity, the solution takes the form 



COS z 



+ zVl-z^, (38) 



(l+^0)3/2 

where the radial coordinate is written in terms of the variable z, which is defined by 



2_ 1 + e 



1 + ?0 

Note that equation (j38p determines the period of a radial orbit for a given starting point ^oi 
specifically, the half-period is given by 

1 



ri/2 = ^ [(1 + Co)'/' cos-i y^TTo + V^o + eo\ ■ (40) 

In the inner limit where <C 1, the solution of equation (j38p reduces to the solution of a constant 
force equation and the time scale of a half period (the time required for on inward or outward 
crossing) simplifies to the form 

ri/2 = • (41) 

For general (non-radial) orbits, the half-periods and turning angles have been studied previously 
for the case of a spherical Hernquist potential (AB05). 
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3.3. Brief Survey of Orbits 

The triaxial system considered herein supports a wide variety of orbits, and displays typical 
dynamical behavior for this class of astrophysical system (BT87) . Previous work shows that as the 
axis ratios increase (i.e., as the ellipsoidal density contours become more elongated), the fraction 
of orbits that are chaotic increases (PMOl). Similarly, as the axis ratios increase, the fraction of 
box orbits increases and the fraction of loop orbits decreases (BT87). Both of these trends are 
illustrated in Figure 1, which shows a series of surfaces of section with increasing axis ratios. For 
all three systems, the maximum value of the triaxial coordinate TTl is taken to 1)6 TTZjxiax — 

0.1; the 

maximum value of the x-coordinate is then given by x^^ax = arriraax (where a is the value of the 
axis weight). The energy measured relative to the minimum value at the center of the potential 
well is nearly constant (with values AE = 0.205, 0.224, and 0.216 for the cases shown). Since we 
are primarily interested in the instability of orbits that are initially confined to one of the principal 
planes, these orbits are all fixed in the x — z plane; similar results hold for orbits started in the other 
two planes. The first plot in Figure [Ta| uses axis weights near unity so that the mass distribution 
is nearly spherical. As expected (BT87), the surface of section for this system is well ordered and 
most of the phase plane corresponds to loop orbits. Note that loop orbits correspond to curves in 
the surface of section where the end points terminate along the positional (x) axis, whereas the box 
orbits have patterns in the plane that continue into the velocity (vertical) axis. As the axis weights 
depart farther from unity (see Fig. [Ta|) . the surface of section plot breaks up into an increasing 
number of islands and contains a lower fraction of loop orbits. 
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Fig. la. — Surfaces of section for triaxial potential with varying axis weights, (a) Surface of section 
for nearly spherical version of the potential with (a, c) = (VTT, \/(L9). (b) Surface of section for 
moderate axis weights (a, c) = (1.25,0.75), so that the axis ratios are 3:4:5. (c) Surface of section 
for elongated potential with (a, c) = (\/3, \/2/2). The letters refer to the orbits shown in Figure [2j 
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Fig. 1 — Continued. 
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Fig. 1 — Continued. 

To illustrate the variety of possible orbits, the third surface of section in Figure [la] includes 
labels for 15 points in the phase plane (indicated by the letters 'a' - 'o'). Each of these points is 
associated with a corresponding orbit, which is shown in the 15 panels of Figure [21 Orbit 'a' is 
located in the central region of the phase plane, i.e., the region that dominates for nearly spherical 
potentials (see Fig. [Taj); the corresponding orbit in Figure [5] is thus a tube orbit. The tube orbits 
become increasingly complicated (see orbits 'b' and 'c') as one moves out from the central region. 
The remaining parts of the surface of section primarily correspond to box orbits of varying degrees 
of complexity. For orbits that lie in the centers of islands in the phase plane, the corresponding 
orbits show symmetry ~ for example, orbit 'i' is a 'fish orbit'. Notice that orbits that lie at the 
centers of islands in the phase plane exhibit this symmetry because they represent "resonant" 
orbits. For example, the fish orbit in Figure [2] corresponds to a 2:3 resonance in the x — z plane. 
Note that the orbits in Figure [2] were constructed for purposes of illustration, so that "tighter" fish 
orbits, closer to exact resonance, can be studied. One should also keep in mind that these orbits 
are confined to the x — z plane; additional types of orbits can be obtained when the orbit is allowed 
to explore the full three dimensions of the potential (BT87, PMOl). The instability of these orbits 
to motion in the perpendicular direction is considered in the following section. 

Figure [3] shows a typical box orbit in this triaxial potential. The axis ratios are chosen to be 
moderately triaxial, but not extremely so. This orbits displays a combination of regularity and 



Fig. 2. — Collection of orbits for triaxial potential with axis parameters (a, c) = (\/3i V^/'^)- The 
panels are labeled by letters 'a' - 'o', where the labels coincide with those on the surface of section 
plot shown in Figure Ic. 
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chaotic behavior. The orbit crosses back and forth across the origin with no well-defined sense 
of rotation. The crossing time scale is of order unity (in dimensionless units), although it varies 
slightly from cycle to cycle. The distance of closest approach to the origin also changes from cycle 
to cycle, but exhibits greater variation. In the discussion of orbital instability that follows, we use 
this "typical" orbit as a benchmark case to study the instability. 



4. INSTABILITY 

A numerical survey of parameter space reveals that many of the orbits that are initially con- 
fined to one of the principal planes are unstable to motions in the perpendicular direction. One 
example of this unstable behavior is shown in Figure [H In this case, the orbit begins in the x — z 
plane, but has a small perturbation in the y-direction (here Ay = 10~*). As shown in Figure [H 
the subsequent evolution of the y coordinate displays four types of behavior: [1] The most striking 
property of the solution is the overall growth of the amplitude, i.e., the envelope shows exponen- 
tial growth. [2] The solution exhibits a well-defined periodicity superimposed on the exponential 
growth. This periodicity is a reflection of the underlying (near) periodicity of the original box orbit 
in the x — z plane and is thus not unexpected. However, notice that the periodicity has multiple 
frequencies. [3] The solution displays chaotic irregularities superimposed on the otherwise smooth 
(regular) function. [4] The growth saturates at a well-defined time scale, and the solution subse- 
quently undergoes long period oscillations (with additional short period oscillations and chaotic 
irregularities superimposed). In the analysis and discussion given below, our goal is to elucidate 
the fundamental mechanism(s) responsible for these four properties. 

For orbits that begin in a principal plane, the equation of motion in the perpendicular direction 
(out of plane) determines whether or not the orbit is stable. For a given plane, we can consider the 
perpendicular coordinate value to be small, at least at the beginning of the evolution. For example, 
if the original orbit lies in the y — z plane, then the x coordinate is small so that |x| <C y, z. In this 
limit, the equation of motion for the non-planar direction takes the form 

— +u;lx = where ^ J , (42) 

Similarly, for the other two principal planes, the equations of motion for the perpendicular coordi- 
nate can be written 

— 2 + Wj^y = where uj = —===== — , (43) 



and 



d^z 9 ,9 4/c , , 

—J + uJlz = where = = ^^^= . (44) 

dt ^JlP'x'^ + a^y"^ + cy x"^ + 
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Fig. 3. — Box orbit initially confined to the x — z plane. The potential is taken to have axes 
{a,b,c) = (\/2, 1, 0.5071). For this orbit, the time average (j) of the specific angular momentum 
approaches zero, so the orbit has no sense of rotation and is thus a box. The orbit crosses near the 
origin in a quasi-periodic manner, which leads to quasi-periodic forcing of the perturbations in the 
perpendicular direction (see text). 
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Fig. 4. — Development of the instability for the box orbit shown in Figure El The original orbit 
lies in the x — z plane with a small perturbation in the perpendicular direction (Ay ~ 10~^). 
The amplitude of the perpendicular coordinate (the y coordinate) thus starts at \y\ ~ 10^^ and 
grows to saturation during ~ 100 time units. Note that we take the absolute value of y{t) and 
that we scale the coordinate by its initial value, so the vertical axis represents the growth factor. 
Notice also that the time evolution of the perpendicular coordinate simultaneously displays nearly 
periodic behavior, exponential growth, and chaotic variations; at late times the exponential growth 
saturates and is replaced by long-period oscillatory behavior. 
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Fig. 5. — Driving frequency as a function of time for the orbit depicted in Figures [3] and 

m This plot represents the effective oscillation frequency of the test particle in the direction 
perpendicular to the starting plane of the orbit (here, the y direction). The three panels show three 
temporal ranges, i.e., 1000 time units (top), 100 time units (middle), and 10 time units (bottom). 
Notice that the spikes in u;^ are extremely narrow, much closer to the delta function limit than the 
cosine limit. Notice also that the forcing function has structure on a wide range of time scales, and 
shows both periodicity and chaotic behavior. 
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Fig. 6. — Histograms showing the distributions of periods (top panel) and forcing strengths (bottom 
panel) for the orbit considered in Figures 2-4. The period is defined to be the time interval between 
successive maxima in the function io'^{t) and the forcing strength q is defined to be the peak values 
of the function (the maxima). Note that the scale for the period, and hence its distribution, is 
narrow, with mean value (Porb) = 0.986, and standard deviation ap = 0.057. In contrast, the 
distribution of forcing strength is shown using logarithmic binning and spans a wide range, with 
mean value (g) ^ 360 and standard deviation ag ^ 560. 
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These equations of motion for the perpendicular coordinate (eqs. |42H44j ) provide both a 
quahtative and quantitative explanation for the instability, although the rest of the orbit, described 
by the other coordinates [e.g., x(t) and z{t)], play a role. (In addition, in order to understand the 
saturation mechansim, we must go beyond the limit of small perpendicular coordinate values — see 
below.) To fix ideas, we present the driving frequency tOy for a box orbit in the x — z plane in Figure 
[5l This function io'^{t) is typical for orbits that exhibit instability in the perpendicular coordinate, 
although the exact form varies with the orbit under consideration. Note that the driving frequency 
shows both periodic and stochastic behavior. Furthermore, this form for u;^ allows for a qualitative 
description for the four properties of the unstable solutions listed above: 

[1] Instability. At first glance, the equations look like those of a simple harmonic oscillator, 
which would naively imply stability. Upon closer inspection, however, the functions io'^{t) are highly 
time variable and nearly periodic. In the limit where the functions uj'^{t) are exactly periodic, the 
equations of motion are a type of Hill equation (Hill 1886; Magnus & Winkler 1966, hereafter MW66; 
Abramowitz &; Stegun 1970, hereafter AS70). Since Hill equations are known to have regimes of 
instability, the fact that we see unstable solutions follows immediately. Notice that in the limit of 
nearly circular orbits in a nearly spherical potential, the equation of motion for the perpendicular 
coordinate becomes Mathieu's equation (§5.3), whose regimes of instability are well studied (AS70, 
Arfken & Weber 2005, hereafter AW05). This Mathieu-equation limit has been studied previously 
in a related context (Binney 1981, Scuflaire 1995). Notice also that the driving frequency function 
uj'^{t) displays extremely narrow spikes, rather than smooth cosine curves, as in Mathieu's equation, 
so the latter does not represent a good approximation for most of the parameter space of interest 
here. We thus consider the opposite limit in which the peaks in the driving frequency function are 
considered as Dirac delta functions (see §5; Appendices B and C). 

[2] Periodicity. The second key feature of the observed unstable solutions is their periodicity. 
This feature also follows directly from the fact that the perturbation equation is of the Hill variety. 
In particular, Floquet's Theorem implies that the Hill equation has solutions of the form 

y{t) = e-*p(t) , (45) 

where p{t) is a periodic function and where the phase factor can provide a growing (unstable) 
solution when a is imaginary (MW66, AS70). These previous mathematical results imply that the 
period of p{t) is either the same as that of the driving frequency w^(t) or twice as great. In this 
context, the periodicity of the driving frequency uj'^{t) is determined by the period of the orbit 
in the (original) orbit plane. This period, in turn, is given by equations (|40l) and (l4T]) to a good 
working approximation. As shown in Figure O and indirectly by Figure O the driving frequency 
has multiple periods; this complication leads to the solutions y{t) displaying multiple periodicities 
as well. 

[3] Chaos. The motion of the perpendicular coordinate is ultimately driven by the (time 
dependent) coordinates of the starting orbit (in its initial plane). The amplitude of the driving 
frequency oscillations uj'^{t) is essentially determined by the distance of closest approach of the orbit 
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to the origin (see eqs. [l2]-[33]) and this distance varies from passage to passage. If the orbit is 
chaotic, as is often the case, these distances of closest approaches wiU be given by some distribution 
of values, and the corresponding driving frequency will have a chaotic element. 

One might worry that departures from exact periodicity would invalidate previous mathemat- 
ical results concerning instability. This issue led us to analyze the quasi-periodic Hill's equation 
with random forcing strengths (Adams & Bloch 2007, hereafter AB07). These new results show 
that the departures from strict periodicity in Hill's equation result in a re-scaling of the parameters, 
but do not otherwise compromise the instability (see Theorem 1 of AB07; see also Appendix D). In 
addition, the stochastic variations in the forcing function lead to two separate contributions to the 
growth rate: an asymptotic growth rate that results from a direct application of Floquet's theorem, 
and an anomalous growth rate that results from matching the solutions from cycle to cycle, where 
each cycle has a different forcing strength (see Theorem 2 of AB07). 

The distributions of periods and forcing strengths thus play an important role in determining 
the dynamics. For the box orbit depicted in Figure [3l and the corresponding forcing function (^"^{t) 
shown in Figure \5\ we show the distributions of periods and forcing strengths in Figure [6j For 
the sake of definiteness, we have defined the period to be the time interval between peaks of the 
function w'^(t). The forcing strengths (q) are defined to be the heights of the peaks. Although 
the distributions of period and forcing strength vary with the orbit under consideration, Figured] 
reveals several important and typical trends: The distribution of periods is quite narrow, so that 
departures from strict periodicity are small. This finding makes the corrections for varying periods 
corrrespondingly small (where this statement is quantified in Appendix D). On the other hand, 
the distribution of forcing strengths q is wide, with the standard deviation larger than the mean. 
Both the mean value (appropriately weighted) and the width of the distribution play a role in 
determining the growth rates for instability. 

In addition to chaotic behavior within a given orbit, the variation of dynamical properties 
varies from orbit to orbit, even those with almost the same starting conditions. The orbits thus 
display sensitive dependence on their initial conditions, one of the defining properties of chaotic 
systems. As one example of this behavior. Figure [7] shows the growth rate of the instability plotted 
as a function of the starting x-coordinate. This plot shows that nearby starting trajectories can 
result in significantly different growth rates. Since the growth of the y-coordinate is not perfectly 
exponential, some ambiguity arises in the definition of the growth rate. The results shown in Figure 
[7] were obtained by finding the maximum value |y|max of the y-coordinate over the first 50 time 
units, and defining 7 = In \ ymax\/50; these growth rates thus represent the average growth rates over 
this time interval. Figure [7] also shows the saturation level of the perpendicular coordinate (in the 
lower panel) plotted as a function of the same starting conditions. Here, the saturation levels were 
computed by finding the maximum displacement over the first 200 time units (thereby allowing the 
instability time to saturate). The behavior of the saturation level is also highly structured, i.e., not 
a smooth function of XQ/xmax- The fine scale structures in both the growth rate and the saturation 
level (when considered as a function of xo/xmax) arise (in part) from the narrow bands of stability 
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that are always present in the Hill's equation that describes the instability (see AS70, Weinstein 
& Keller 1987, AB07, and §5 below). Further, we find that the saturation levels correlate with the 
growth rates; in particular, low or vanishing growth rates correspond to low saturation levels, as 
expected. Figure [7| shows additional systematic behavior, in that the ranges of starting conditions 
that lead to large growth rates and large saturation levels are relatively well-defined. The larger 
structures in the growth rate functions (and saturation levels) often coincide with the structure of 
the phase-plane diagram for this potential. [4] Saturation. The strength of the instability, as well 
as its existence, depends on the amplitude of the driving frequency u;^(t). The equations of motion 
found above are derived in the limit where the perpendicular coordinate is small compared to the 
coordinates in the (initial) orbit plane. As the perpendicular coordinate grows, this approximation 
becomes invalid, and the relevant equation of motion takes on a modified form. For the cases 
of interest here, the driving frequencies have the basic form w^(t) ~ 1/^,, where is a modified 
"radial coordinate" where the (x, y, z) coordinates are weighted differently, by the geometric factors 
(a, 5, c). As the perpendicular coordinate (e.g,. y) grows, its contribution to and hence the driving 
function becomes significant, with the net effect being to reduce the forcing amplitude. When the 
amplitude decreases enough, the motion is no longer unstable in the perpendicular coordinate, and 
the instability saturates. 

As a particular example of the behavior described above, we consider the equation of motion 
for the y coordinate for an orbit initially confined to the x — z plane in the limit of a nearly spherical 
potential (eq. [22]), i-e., 

d?y , 1 r, , K^'^ - z^) 



dt^ ^ 2i 



1 + 



y = 0. (46) 



4^2 

When the y coordinate is small (compared to x and z), the size of the driving term is determined 
by 1/C ~ {^'^ + -z^)"^''^- The forcing term is large for small values of displacement ^. When the y 
coordinate grows, however, its contribution to ^ becomes significant (i.e., ^ is not necessarily small 
even when |2:| ^ 1), and the instability saturates. 

To see this trend explicitly, let = {d^ + t^y^'^, where this form assumes that the orbit 
has a distance of closest approach d and the radius of the orbit is given by a constant force law 
(appropriate for radial orbits - see §3.2). For a nearly linear orbit with closest approach d, the 
radial distance is given by ^'^ = + , where s is the distance along the orbit measured from the 
point of closest approach. Since s (xt"^ for a constant force law, one obtains the relation = d'^ + t^ 
in appropriate dimensionless units. The effective amplitude of the driving frequency is then given 
by the integral 

dt ^ 8r2(5/4) ^ 371 

(d2 + t4)l/2 ^ Vd' ^ ^ 

At the start of the evolution, when the motion is confined to a plane, the distance of closest approach 
to the origin d = do is given by that of the starting orbit. For a box orbit, do can be quite small, 
so that q is large. A general property of Hill's equations is that larger q implies greater instability 
(although islands of stability will remain - see MW66, AS70, and §5 below). Once the instability 
has developed, the distance of closet approach must include the displacement in the perpendicular 
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Fig. 7. — The upper panel shows the growth rate of instabihty as a function of starting coordinate. 
For all of the cases considered here, the energy is taken to be = 0.1, and the axis variables b = 
1 and c = ^/2/2. Solid curve shows the growth rate 7 for varying starting locations for the axis 
parameter a = -^3/2; the dashed and dotted curves show the results for a = \/3 and a = \/6, 
respectively. The bottom panel shows the corresponding saturation level, i.e., the maximum value 
of the y coordinate expressed as a fraction of the typical x coordinate. The three curves correspond 
to the three values of a in the top panel. 
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direction, so that d — > {y) (where the angular brackets represent an appropriate average of the 
perpendicular coordinate). The value of q at later times is thus reduced by a factor (do/(y))^^^- 
Notice that saturation of the instability requires that (y) ^> do, which is much less restrictive than 
the requirement (y) ~ dmax- In other words, the instability can saturate when the perpendicular 
displacement, measured here by (y), becomes sufficiently larger than the (initial) distance of closest 
approach, but need not be comparable to the full orbit size. In many cases, however, the saturation 
levels are large enough that (y) is comparable to the original orbital extent in the x — z plane (see 

Fig. ED. 
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Fig. 8 — Continued. 
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Fig. 8a. — Instability as a function of axis weights. Nearly spherical potentials correspond to the 
upper right corner of each panel. For each point in the grid, orbits are started in the x — z plane, 
with fixed energy, and with a fixed starting value of the x coordinate. The symbols denote the type 
of resulting orbits, where open squares are stable box orbits, open triangles are stable loop orbits, 
stars are unstable box orbits, and crosses are unstable loop orbits, (a) Result for xo/xmax = 0.10. 
(b) Result for xo/x^ax — 0.50. (c) Result for xq/x^ux — 

0.90. 
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Fig. 8 — Continued. 

Before leaving this section we consider the degree of stabihty versus instabihty as a function 
of the shape of the triaxial potential, where shape is codified by the axis weights (a,6, c). For a 
given energy and set of axis weights, we start orbits with a given value of x^/xmax as described in 
§3. We then determine if the orbit is a loop orbit or a box orbit, and whether the orbit is stable 
or not. To determine the orbit class, we use a procedure similar to that developed by Fulton & 
Barnes (2001). For a given orbit, we fourier analyze both time series x{t) and z{t) and obtain the 
dominant frequencies. If the dominant frequencies are in 1-1 correspondence, the orbit is labeled 
as a loop (where we allow for a 5% discrepancy in the frequency ratio). To determine if the orbit 
is unstable, we monitor the maximum displacement \y\M in the perpendicular direction. If the 
perpendicular coordinate grows above a given threshold, then the orbit is considered unstable. 
Since orbits generally either grow and saturate at large values of \y\M, or remain stable with small 
values of y ~ 10~*, the results are insensitive to the threshold value. For the sake of definiteness we 
use a threshold \y\M/fM ^ 0.01 (where vm is the maximum two-dimensional radius in the original 
X — z plane of the orbit). The results are shown in Figure [Ha] for three different values of the 
starting x coordinate. Several trends are evident. In all cases, a significant fraction of the plane 
shows unstable orbits, so that the instability considered herein is indeed robust. As the starting 
value xa/xmax is increased, the orbits become more nearly "radial" (following the x axis) and the 
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degree of instability decreases somewhat (compare the panels of Fig. [8a|) : further, the regimes of 
stability (usually) correspond to extreme axis weights (small c and large a). Finally, we note that 
the regimes of stability and instability are intermixed, as expected in highly chaotic systems. 

5. MODEL EQUATIONS 

In this section we develop and analyze a class of model equations that captures the essential 
physics of the orbit instability found above. As shown in the previous section, when the motion 
is initially confined to one of the principal planes, the equation of motion for the amplitude of the 
perpendicular coordinate takes the form 

5 + [A + Q(t)]y = 0, (48) 

where Q{t) is nearly periodic. Here we use the variable y as the perpendicular coordinate, but 
the same general form holds for the perturbation direction falling along any of the axes. Note 
that we have separated out part of the forcing function as a baseline frequency ^/X, although this 
decomposition is not unique. Notice also that the function Q{t) can be considered as a function 
of the form Q[x{t), z{t)], where x{t) and z{t) describe the orbital motion in the original plane, so 
that Q is given by equations (fT8H20|) in general and by equations (P2ti3^ in the limit of small y 
(as assumed here). 

In most applications, Q{t) will be nearly periodic because its time variation is given by an orbit 
solution, although the crossing time of the orbit can vary from cycle to cycle. In this context, in the 
Hernquist potential, the orbital periods depend on energy, but are nearly independent of angular 
momentum (see Fig. 3 of AB05), especially in the inner limit as considered here, so that the cyclic 
forcings have nearly the same time intervals. This expectation is borne out in numerical evaluations 
of the orbits, as shown in the distributions of Figure [H On the other hand, the amplitude (the 
forcing strength) varies substantially from cycle to cycle. If the function Q{t) is exactly periodic, 
then the relevant equation of motion has the form of Hill's equation (Hill 1886), and we can draw 
on a number of well-known results (MW66, AS70, Arnold 1978). In this section we consider the 
periodic version of Hill's equation as a first approximation, and then show how to generalize the 
results to cases where the forcing strength varies from cycle to cycle. This section is augmented 
by Appendices B - D, which include some of the mathematical details, and by Appendix E, which 
presents a heuristic review of Hill's equation. We consider a more rigorous mathematical treatment 
of the stochastic problem in a separate publication (AB07). 

5.1. The Delta Function Limit 

As a first approximation, we consider the forcing potential to be sufficiently sharp that we 
can model Q{t) with a Dirac delta function. In physical terms, this approximation means that the 
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orbit passes close to the origin (so that the spike in Q{t) is large) and that the orbit spends little 
time near the origin (so the spike is narrow). Both of these conditions are met for box orbits in 
triaxial potentials. In this limit, the equation of motion for the perpendicular coordinate takes the 
particular form 

^ + [X + q5{[t]-7:/2)]y = 0, (49) 

where q measures the strength of the spike and where 6 is the Dirac delta function. Here we 
have scaled the time variable so the period of one cycle is vr. Note that the argument of the delta 
function is written in terms of [t], which corresponds to the time variable mod-vr, so that the forcing 
potential is vr-periodic. Notice also that the actual physical problem has a forcing strength q that 
varies from cycle to cycle. As a result, we should think of the "constant" q as a stochastic variable 
that takes on a new value every cycle < [t] < vr. 

In the delta function limit, the solution to Hill's equation is thus specified by two parameters: 
the frequency parameter A and the forcing strength q. The solutions to equation (j49p are constructed 
in Appendix B along with the determination of the growth rates for instability. Figure [9] shows 
the plane of possible parameter space for Hill's equation in this limit, with the unstable regions 
shaded. Notice that a large fraction of the plane is unstable. This result is in keeping with our 
previous finding that the equation of motion for the instability takes the form of a Hill's equation 
with highly spiked forcing functions and that a large fraction of the orbits are unstable (§4). 

An important feature of these solutions is that in the limit of large forcing strength q ^ 1, the 
bands of stability become extremely narrow — but they never vanish entirely. Physical orbits are 
often found in this limit, as illustrated by Figures[5]and[6l This property is generic to Hill's equations 
in the unstable (here, large q) limit (Weinstein &; Keller 1987). For the case of delta function forcing 
terms, for example, the widths of the stable zones are given by (AA) = (8n^/7r)(l/g), where n labels 
the order of the stability strip (AB07) . The presence of these narrow zones implies that the issue of 
stability /instability can depend sensitively on the orbit parameters, i.e., nearby orbits can display 
significantly different behavior. This trend is seen in the numerical results for the growth rates, 
as depicted in Figure El where the function 7(xo/xmax) displays a great deal of structure, so that 
orbits with nearly identical starting conditions can have significantly different growth rates. 

With the solution to the model equation in hand, as illustrated by Figure O it is useful to 
reconnect with the original orbit problem. In terms of physical orbits, the period of the (periodic) 
forcing potential is determined by the crossing time tcross of the orbit. Here we have defined the 
dimensionless time variable t so that the period is vr; the physical time variable tphys is thus related 
to the dimensionless time appearing in equation (jl9]l via tphys = t{tcross/'^)- Further, the period of 
a half orbit for the spherical limit is given (almost) analytically in AB05, 

Ti/2 = (47rGpo)~^/^e"^/^ {cos + V^Vl^} A , (50) 

where e = |i?/^'o| is the dimensionless energy and where the correction factor A lies in the range 
1 < ^ < 1.05 (see Fig. 3 of AB05). The strength q of the forcing potential is determined by the 
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distance of closest approach of the orbit to the origin {q ~ l/d). In the spherical limit, the distance 
of closest approach is given by the inner turning point of the orbit and is a simple function of energy 
and angular momentum (AB05). Here, for loop orbits, the orbit has an analogous inner turning 
point and the spherical results can be used as a good starting approximation. For box orbits, the 
test particle wanders arbitrarily close to the origin, so that the distance of closest approach will 
often be small (and hence the parameter q will be large); in general, q can vary from cycle to cycle. 
The problem of Hill's equation with random variations in forcing strength is addressed in Appendix 
C for the delta function limit (see AB07 for a more general treatment). 



5.2. Beyond the Delta Function Approximation 

The physical orbits do not produce perfect delta functions for the forcing potential Q{t). 
Instead, the spikes have a small but finite width. In order to understand how this feature affects 
the growth rates (in particular, the regions of stability and instability), we model the forcing 
function through an equation of the form 

+ [A + g„sin-t]e = O, (51) 

where n is a (generally large) even integer. For the class of orbits considered in the previous section, 
we have fit functions of the form sin" t to the functions iv'^{t) resulting from numerical integration of 
the equations of motion. The values of n required to fit these results fall in the range n = 10 — 100, 
with a "typical" best-fit value n ~ 64. The full width at half maximum (FWHM) of the spike is 
given hy w ^ 2([2 In 2]/n)^/^ in the limit of large n where the spikes have narrow widths. For n = 
64, for example, we obtain w ~ 0.294. 

To compare model equations of this form for different values of n, or to compare with the delta 
function limit, one needs to define the effective forcing strength for a given value of qn- In the delta 
function limit, this term has the form q6{[t] — 7r/2). Since the delta function integrates to unity 
over one cycle, we can define the effective strength geff for a given n through the relation 

qeS = qn sm dt = qnTT . (52) 

Jo nil 



Figure [TO] shows the plane of parameters {X,q(,ff) for a Hill's equation of the form given by 
equation (j5ip with n = 64. The q values have been scaled according to the normalization of equation 
()52p so that this plane can be directly compared with Figure [H After including this normalization, 
the two planes are more similar than different, so that we expect the delta function model of 
the previous subsection to provide a good description of the instability mechanism. Nonetheless, 
Figures [91 and [TOl show slight differences, with the "tongues of instability" being somewhat narrower 
for the case with uj'^ spikes of finite width. 
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Fig. 9. — Regions of instability for Hill's equation in the delta function limit. The shaded regions 
show the values of {X,q) that correspond to exponentially growing (unstable) solutions, which 
represent unstable growth of the perpendicular coordinate for orbits in our triaxial potential that 
are initially confined to one of the principal planes (see also AB07). 
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Fig. 10. — Regions of instability for Hill's equation for periodic spikes with nonzero width. The 
width is given by sin"f where n = 64. The shaded regions show the values of {X,qeff) that 
correspond to exponentially growing (unstable) solutions (compare with Fig. E]). Note that the 
values of q^f / have been normalized according to equation (|52|) . 
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5.3. The Nearly Spherical (Mathieu Equation) Limit 

We note that previous treatments of similar orbit instabilities (e.g., Binney 1981, Scuflaire 
1995) often find that the equation of motion for the perpendicular coordinate (whose motion can 
be unstable) takes the form of Mathieu's equation 

+ (a + qcost)y = . (53) 

Although this equation has qualitatively similar instabilities to those considered here (compare 
the results presented in AS70 with Figures [9] and [T0\ see also Binney 1981), the forcing function 
Lo'^ ~ cos t is extremely different from that shown in Figure [5l It is thus useful to find the limiting 
regime of parameter space for which the instability takes the form of Mathieu's equation, rather 
than the more extreme, nearly delta function limit considered above. 

In qualitative terms, in the limit of a nearly spherical system (§2.5), the equation of motion 
for the perturbation takes the simpler form 

where ^'^ = x'^ -\- y'^ + and where we have assumed motion initially confined to the x — z plane. 
If we also consider the initial orbit to be nearly circular, then we can define an effective radial 
coordinate of the orbit = + . The forcing term thus takes the form 

2 _ 2 _ 1 _ 1 
^ ~ l'^ rit) ^ ro + (Ar)cos$7t ' ^ ^ 

where ro plays the role of the semimajor axis of the orbit, the turning points of the orbit fall at 
rgib Ar, and O is the mean motion of the orbit (since the orbit is not an ellipse, these quantities have 
slightly different meaning than in the standard Kepler problem - see AB05 for further discussion). 
With one further approximation, assuming Ar <^ ro, the equation of motion for the instability 
takes the form 

d^y 2 r Ar ^ 1 

-4 + — 1 cos m y = , (56) 

at^ ro L ro J 

which has the form of Mathieu's equation. Although this result was derived for a nearly spher- 
ical potential, the loop orbits in triaxial potentials also result in Mathieu-like equations for the 
instability of the perpendicular coordinate (although the fraction of loop orbits tends to decrease 
with increasing levels of triaxiality) . Note that Mathieu's equation can also be transformed into a 
discrete map using the same procedure as described in the previous subsections. 



6. CONCLUSION 



6.1. Summary of Dynamical Results 

Analytic results: One set of results from this paper is the development of a collection of 
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analytic expressions that describe potentials, forces, and orbits for the triaxial generalization of the 
Hernquist/NFW profile in the inner limit. In particular, we have found a completely analytic form 
for the potential (eqs. [I0]-[T7]), the force terms (eqs. [18] -[20]; see also PMOl), the axis ratios 
of the potential (eqs. [Hi-El]), and the potential in the limit of nearly spherical density profiles 
(eqs. [27] -[28]). We have also found the radial solutions in the spherical limit for the full Hernquist 
potential (§3.3), i.e., valid over all radii < ^ < oo, not just the inner limit ^ 1. Finally, for the 
axisymmetric version of the NEW profile, we have found analytic expressions for the potential and 
the force terms that are valid over the entire radial range of the halo (see Appendix A). 

Another contribution of this paper is an analytic formulation of the instability problem. For 
orbits initially confined to one of the principal planes, the equation of motion for the perpendicular 
component takes a relatively simple form (eqs. [12] -[44]). Further, the instability can be described 
in terms of simple model equations (§5, Appendices B - D) that allow for an analytic description 
of the instability. Since these equations of motion have the form of Hill's equation, we can utilize 
a large collection of known results (AS70, AW05, MW66). In this context, however, the relevant 
Hill's equations have random forcing terms (forcing strengths that vary from cycle to cycle), which 
requires the development of new mathematical results. We have addressed this complication in a 
separate companion paper (AB07). 

In addition to their usefulness in providing physical understanding of orbits and their instabil- 
ities, these analytic results are also useful for computational purposes. Although all of these results 
can also be determined numerically, the regime of parameter space is enormous, so that analytic 
results provide a substantial savings in the required computational effort. Analytic results also aid 
in our understanding of the underlying mechanism for the instability, as outlined below. 

Overview of the instability: The main focus of this paper is to demonstrate that for orbits 
initially confined to any of the principal planes, the motion in the perpendicular direction is unstable 
to growth; our related objective is to understand the underlying mechanism for the instability. 
As discussed in §4, the instability displays four key features: quasi-periodic oscillatory behavior, 
exponential growth, superimposed chaotic variations, and eventual saturation. These four features 
can be understood as follows: 

For a given energy, an orbit has a well-defined crossing time, which naturally gives rise to 
quasi-periodic behavior. Furthermore, this property provides an effective periodic forcing potential 
for the motion in the direction perpendicular to the original orbital plane. The resulting dynamics 
can thus be described with a Mathieu-Hill-like equation, which is well-known to have unstable 
(growing) solutions (MW66, AS70, AW05). As a result, for a wide range of parameters, the motion 
in the perpendicular direction will be unstable. In the case of a "pure" Hill's equation, with 
precisely periodic forcing, the growing solutions have exponential growth superimposed on periodic 
oscillations, as seen in the instability under study. The period of the growing solutions is given by 
the period of the forcing function, which is determined here by the crossing time of the orbit in the 
original plane. Thus, the first two properties of the instability (periodic behavior and exponential 
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growth) follow directly from the form of the equation of motion for the perpendicular coordinate 
(eqs. [12]~[11]). 

Chaotic variations represent the next property of the instability and they arise due to the 
triaxial nature of the potential (which arises from the triaxial density profile). As is well known 
(BT87), triaxial potentials allow for box orbits, and for chaotic behavior, and both of these effects 
are increasingly common as the axis ratios of the triaxial potential become more extreme. For 
sufficiently triaxial potentials, the crossing time varies from cycle to cycle and hence the frequency 
of the forcing of the instability varies somewhat. However, a much more significant effect is that the 
strength of the forcing varies from cycle to cycle. In this case, the forcing strength is determined 
by the distance of closest approach of the orbit to the origin in the original plane, where this 
distance is weighted by the axis weights (a, 6, c) of the triaxial potential. For example, box orbits 
are quasi-periodic and wander arbitrarily close to the origin. Thus, as the orbit varies chaotically 
in its original plane, the forcing terms for the instability vary as well, and these chaotic variations 
appear in the growing (unstable) solution for the perpendicular coordinate. 

The final property of the unstable solutions is their eventual saturation. This behavior is 
expected for two reasons: In physical terms, the full three dimensional orbit must conserve energy, 
so there is a limit to the amplitude of oscillations in the perpendicular (unstable) direction, i.e., 
saturation is inevitable. In mathematical terms, the forcing strength depends (inversely) on the 
distance of closet approach to the origin; for orbits with sufficiently large values of the perpendicular 
coordinate, the forcing strength becomes small and the instability stalls. 

Model equations: We have analyzed the instability of Hill's equation in the delta function limit, 
where the forcing terms are localized in time. In this limit, the equation of motion can be solved 
analytically for a given cycle of the motion in the original plane. As a result, the growth rates for 
a pure Hill's equation can also be found analytically (§5 and Appendix B). For the case of forcing 
strengths that vary from cycle to cycle, we have found estimates for the growth rates in terms of 
the expectation value of the growth rates for individual cycles (Appendix C). In the limit of highly 
unstable systems, this expectation value - the asymptotic growth rate - provides a good working 
estimate. However, the true growth rate contains an additional component (the anomalous growth 
rate) that arises from the matching conditions at the cycle boundaries. In a separate paper we show 
that this correction term can be found analytically in terms of an expectation value (see AB07), 
and we present general theorems regarding the instability and growth rates for Hill's equation with 
random forcing terms. 

The y direction is most unstable: Although orbits around all three axes are unstable, the orbits 
around the intermediate (y) axis are more likely to be unstable: The instability we are exploring 
here is often strongest when the motion in the original plane arises from box orbits, which are more 
likely to take place when the axis ratios of the two axes in the original orbit plane are more extreme 
(BT87). The "y instability" takes place when the original orbit lies in the x — z plane, and the 
axis ratio a/c is the largest by definition. Therefore, the y instability is the most likely to occur; 
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alternatively, the growth rate for the y instabihty is hkely to be larger than that for the other two 
cases. One should compare with situation with the classic result from rigid-body dynamics: As 
is well-known0 rotation about the shortest or longest axis of an asymmetric rigid body produces 
stable motion, whereas rotation about the intermediate axis is unstable. 

6.2. Astrophysical Implications 

Although this paper focuses on the description of the triaxial potential, its orbits, and the 
corresponding orbital instabilities, these results can be applied to a number of astrophysical systems 
on a variety of size scales. On the largest scales, this work applies to the inner regions of dark matter 
halos - both for galaxies and galaxy clusters. On smaller scales, these results are applicable to the 
dynamics of tidal streams, galactic bulges and warps, and young embedded star clusters. All of 
these systems (or the inner portions of them) are described by the triaxial potential developed 
herein, and can be potentially affected by this orbit instability. These applications are briefly 
outlined below, but we emphasize that they should be explored in greater depth. 

For completeness, we note that in real astrophysical systems, the gravitational potential will 
never be completely smooth as assumed here. Instead, stars and other inhomogeneities provide a 
nonzero level of graininess to the potential, and will thereby contribute to the degree of chaos in the 
constituent orbits. In terms of the instability considered here, these irregularities guarantee that 
orbits cannot be exactly confined to any given plane; instead, these effects act to set the smallest 
possible amplitude for perturbations in the direction perpendicular to the starting plane. We also 
note that the instability considered here strictly applies only to orbits that start in one of the 
principal planes, and such orbits represent only a small portion of phase space. 

Dark Matter Halos: In dark matter halos, both the orbits of dark matter particles themselves 
and the orbits of sub-halos benefit from an analytic description, which can aid in understand- 
ing how dark matter halos achieve the (nearly) universal form observed in numerical simulations 
(starting with NFW). Another finding from numerical simulations is that the velocity distributions 
are generally radial in the outer parts of the halo and more isotropic in the inner regions. This 
isotropization cannot take place through two-body relaxation processes, as the time scale is much 
too long. The orbit instability considered in this paper can provide at least a partial explanation: 
The inner regions tend to have triaxial forms with density profiles p ~ 1/m in the inner limit, 
i.e., the form of the density and potential considered in this paper. Since nearly radial orbits are 
subject to the orbit instability considered herein, the inner halo cannot maintain purely radial ve- 
locity distributions. Further, both the time scale (several crossing times) and the saturation level 
(the amplitudes of the perpendicular coordinate oscillations become roughly comparable to those 
in the original orbit plane) suggest that the velocity distribution can become sufficiently isotropic 



This calculation is generally credited to Euler (1749) in classical mechanics texts (e.g., Marion 1970). 
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to explain the results of numerical simulations. The velocity distributions of dark matter particles, 
as well as their enhancements in tidal streams (see below), can affect direct detection strategies 
(e.g., Evans et al. 2000). 

Tidal Streams: Since tidal streams must orbit in galactic potentials with some triaxial aspect, 
the results of this paper also inform their dynamical properties. The coherence of tidal streams 
relies on the orbits of the satellite galaxies and the stars that are stripped away from them to 
remain coherent for a number of orbits. In contrast, the instability studied in this paper can act 
to disrupt tidal streams. This disruption could have important implications for direct searches for 
dark matter, through particle detection experiments, where enhancement due to tidal streams has 
been invoked as a way to help discriminate the signal from the background (Evans et al. 2000). 
Turning the problem around, one can also use the observed properties of tidal streams as a means 
to constrain the triaxiality of their galactic potentials (Ibata et al. 2001). This instability will also 
play a role in the assimilation of merging substructures during the galaxy formation process (e.g., 
Helmi et al. 1999). In all of these applications, however, the number of orbits of the satellites is 
relatively small (several), and it could be hard to distinguish between precession of an ordinary 
orbit in a triaxial potential and this instability. Further exploration of this issues is necessary. 

Galactic Bulges and Disks: The orbits and instabilities studied in this paper can also act on 

the scale of galactic bulges, which can be described by a triaxial generalization of the Hernquist 
potential. In this context, an important issue is how often stars wander close to the potential 
center, where (usually) a supermassive black hole resides. The orbit instability studied herein acts 
to change nearly radial orbits into orbits with more isotropic velocity ellipsoids, and thereby changes 
the rate at which stars interact with the central black hole. We note that near the galactic center, 
the gravitational potential of the black hole itself must be taken into account. Some work along 
these lines has been done (Gerhard & Binncy 1985; PMOl, Poon & Merritt 2002), which shows 
that triaxiality leads to more box orbits and greater rates of stellar accretion by a central black 
hole; on the other hand, the orbit instability studied herein acts to increase the distance of closest 
approach and thereby suppresses accretion of stars. In any case, this issue provides an interesting 
problem for future study. 

A related issue is that sufficiently distorted (non-axisymmetric) disk galaxies can be subject 
to this type of orbit instability, which acts to populate regions out of the original disk plane. This 
issue was explored earlier by Binney (1981) using the triaxial version of the logarithmic potential; 
this work can be used to extend such previous analyses by providing a different potential and 
considering forcing frequencies uj'^{t) that have stochastic amplitude variations from cycle to cycle 
and are highly nonlinear (e.g., closer to the delta function limit than the cosine forcing functions 
that arise for nearly spherical potentials - see §5.3). Both gas and stars can be subject to this 
instability in galactic disks, leading to galactic warps and other interesting structure (e.g., Sparke 
1995, 2002). 

Embedded Star Clusters: The dynamics of young embedded star clusters can also be described 
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using the inner limit of the Hernquist potential (Larson 1985, Jijina et al. 1999, Adams et al. 
2006). These clusters are expected to have more extreme axis ratios than the larger scale systems 
considered above. Observations of embedded clusters (Lada & Lada 2003 and references therein) 
show that the youngest systems are highly irregular, displaying large departures from spherical 
symmetry. Somewhat older systems (a few Myr) are significantly more spherical, so that the process 
of isotropization takes place rapidly. In these systems, orbits can be altered through star-star 
scattering events (where the stars generally have accompanying disks and/or planets) and through 
the instability studied in this paper. If these young clusters form out of initial gas configurations 
that are highly fiattened, this instability will act to populate the regions perpendicular to the initial 
cloud plane. The resulting stellar clusters can thereby become more nearly spherical, even in the 
absence of stellar scattering interactions. For example, suppose that the initial cluster is flattened 
with a 10:1 aspect ratio (c = 0.1). For the orbit instability of this paper to round out the cluster, 
the growth rate 7 must be large enough to compete with stellar scattering. The number Nc of 
crossing times required for the instability to amplify the motion in the perpendicular direction is 
given by Nc ~ (hi 10)/7te- For a purely stellar system, the ratio of the dynamical relaxation time 
to the crossing time (BT87) is given by Nq ~ iV^/(101niV^); for an embedded cluster dominated 
by its gas content, however, the relaxation time is longer so that Nc ~ (Ar^/e^)/[101n(Ar^/e)], 
where e is the star formation efficiency within the cluster (Adams & Myers 2001). These relations 
imply that orbit instability competes with dynamical relaxation when > lOe^ In 101n(A'*/e)/Ar^. 
For typical values (A^^ = 300 and e = 1/3), this constraint becomes > 0.058, a condition is 
that is often met. As a result, orbit instabilities can dominate stellar scattering as a means to 
isotropize young embedded star clusters. In addition, the triaxial nature of the potential in these 
clusters allows young star/disk systems to execute box orbits, which bring the disks closer to the 
massive stars at the cluster center. Triaxial potentials thus allow for greater radiation exposure for 
star/disk systems, compared with spherical clusters, or at least a wider distribution of radiative 
flux experienced by individual solar systems. 



6.3. Discussion and Future Work 

In addition to the applications outlined above, another issue is to consider how this orbit insta- 
bility affects the construction of self-consistent galaxy models. In systems where the orbiting bodies 
provide the source of the potential (e.g., stars in a galactic setting) one must find self-consistent 
models where the orbiting stars provide the gravity for the triaxial potential (Schwarzschild 1979). 
This issue is of particular importance for the intermediate scale of galaxies and galactic bulges, and 
can be complicated by the presence of supermassive black holes (e.g., Poon & Merritt 2002). On 
the larger scale of dark matter halos, the dark matter particles provide the source (mass) for the 
gravitational potential. For these systems, however, N-body simulations of large scale structure 
already self-consistently take into account the orbits of the particles (subject to finite numerical 
resolution) and are the motivation for both the particular forms for the density profiles (NFW, 
Hernquist) and their degree of triaxiality. Finally, on the smaller size scale of embedded stellar 
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clusters, the gravitational potential is dominated by gas, which is subjected to additional forces 
(e.g., magnetic fields) and does not execute the same orbits. As a result, in these systems the 
background potential can be taken as fixed, at least for purposes of studying this orbit instability. 

The analysis of this paper is limited to the case of density profiles that have the form p ~ 1/r 
in the spherical limit, and this form is applicable to only the inner regions of the extended mass 
distributions in question. Future work should consider both the generalization to other inner 
limiting forms (e.g., power-law density profiles with other indices, p ~ l/m"') and extension of the 
potential out to all radii. Toward this latter goal, we present the potential and force terms for the 
axisymmetric version of the NFW profile in Appendix A. In addition to the study of other density 
profiles, the available parameter space for the potentials considered in this paper is enormous, and 
should be explored in greater depth. For example, the case of resonant orbits, including their growth 
rates for instability and saturation levels, would be particularly interesting to consider (anonymous 
referee, private communication). 

Finally, another direction for future work concerns the mathematics of the instability mech- 
anism. This paper presents a class of model equations that describe the instability (§5), and our 
companion paper proves basic Theorems regarding stochastic Hill's equations ( AB07) , with a focus 
on the highly unstable limit. Future lines of inquiry include finding analytic solutions for the case 
of forcing terms with finite width, the relationship between lyapunov exponents for chaotic orbits 
and the growth rates for unstable orbits, better treatment of the saturation mechanism including 
predictions for the saturation levels, and finding (and solving) analogous model equations for gener- 
alized potentials. These mathematical issues, in conjunction with the aforementioned astrophysical 
applications, provide a rich collection of dynamical problems for further study. 
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Appendix A: Potential and Forces for Axisymmetric NFW Profile 

In this Appendix, we present the potential and force terms for the complete axisymmetric 
NFW profile. The results are thus valid over the entire range of radii, not just in the inner limit, 
but they are limited to axial symmetry. 

Evaluating the potential integral for the oblate case, c < 1, yields the result, 

{^/K^+c){l + ^l-q/a){^T^- ^l-q/a) 

(V^ - c)(i - 0^^)(Vr^ + 0^^) 



Va — 7 In 
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+ 2^7 - b 



tan-^ {q/b - 1)"^/^ + tan"^( 
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Note that we have defined a number of new quantities, i.e., 



7 = 1 - 



and 



a,b = - 



1 r 



7 + C'± V(7 + C')'-47i?2 



(A2) 



Keep in mind that a and 6 are the roots defined by equation (A2) and are not the axis weights as 
defined in the main text. For the corresponding prolate case, c > 1, the potential takes the form 
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— — 7 In 
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- c)(l - Vl^^)(vT^ + Vl^^) 

Note that the expressions for the prolate and oblate cases can be put in similar form using the 
identity tan~^x = (i/2)ln|(z + — x)\. Depending on the context, either form can be more 
computationally convenient. 

The forces are found by taking the negative gradient of the potential, F = — V$. Note that the 

potential is taken to be positive in our convention, so that the positive gradient yields the forces. 
For simplicity, we will present the result for the y component of the force; the other forces can be 
found by simple substitution, as outlined below. The force for the oblate case is given by 
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The force for the prolate case is given by 
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The forces in the x and z directions are obtained by replacing the partial derivatives da/dy 
and db/dy with da/dx and db/dx for the x direction, and da/dz and db/dz for the z direction. 
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The above expressions have a number of new variables defined to make them more manageable. 
The partial derivatives of the roots a and b are given by 
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We have also defined the quantities 
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Appendix B: Solutions for the Delta Function Limit 

This Appendix constructs the solutions for the instability equation in the delta function limit 
(see eq. [39]) and finds the growth rate as a function of the input parameters. In this limit, 
the equation of motion has two linearly independent solutions yi{t) and y2{t), which are defined 
through their initial conditions 



yi(0) = l, ^(0) = 0, 



and y2(0) = 0, ^(0) = 1. 

at 



The first solution yi has the generic form 

yi{t) = cos V\t for < t < it/2 , 

and 



yi 



{t) = Acos\/^t + Bsiny/Xt for 7r/2<t<7r, 



(Bl) 
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where A and B are constants that are determined by matching the solutions across the delta function 
at t = 7r/2, including the usual jump condition in the time derivative. We define 9 = \/A7r/2 and 
find 

.2 Q 



A = l + {q/\/X)sm0cos9 and B = -{q/VX) cos" 
Similarly, the second solution y2 has the form 

y2(i) = sin VXt for <t < 7r/2 , 

and 



2/2 



(t) =Ccos^/Xt + DsmVXt for 7r/2<t<7r, 



where 



C = (g/A) sin^ 6* and D = ^ - (g/A) sin6'cos 6* . 
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Using the results derived thus far, we can find the criterion for instability and the growth rate 
for unstable solutions. The growth factor fc per cycle is given by the solution to the characteristic 
equation (see MW66) and can be written as 



fc 



A + (A2 - 4)V2 



where the discriminant A = A(A, q) is defined by 



A = yi(7r) + — (tt). 



(B8) 



(B9) 



It then follows from Floquet's Theorem (see MW66, AS70) that |A| > 2 is a sufficient condition for 
instability. In addition, periodic solutions exist when |A| = 2. At the boundary between stability 
and instability, the parameter A = 2. Solutions with A = 2 are usually unstable, but not always, 
so that further analysis is necessary for this case (MW66). Since the forcing potential is symmetric, 
we know that j/i(7r) = ^2(71"), from Theorem 1.1 of MW66. If we define H = |yi(7r)| = |y2(7r)|, the 
resulting criterion for instability reduces to the form 



H = 



2x/A 



sin(\/A7r) — cos(\/A7r) 



> 1, 



and the corresponding growth rate 7 is given by 

7 = - ln[H + ^/IP^] . 
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Appendix C: Random Variations in Forcing Strength 



The analysis presented in §5 considers the forcing function Q{t) to be perfectly periodic. 
However, the physical orbits found in a triaxial potential often have a random element so that the 
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amplitude q of the function Q{t) varies from cycle to cycle. This Appendix addresses this random 
variation by considering one cycle at a time. Specifically, we consider each period from t = to 
t = vr as a separate cycle, and then consider the effects of successive cycles with varying values of 
forcing strength q. 

During any given cycle (for any value of forcing strength qj), the solution can be written as a 
linear combination of yi and y2- Consider two successive cycles. The first cycle has forcing strength 
qa and solution 

fait) = aaVlait) + Pay2a{t) , (CI) 

where the solutions yia{t) and y2a{t) correspond to those found in the previous subsection when 
evaluated using the value qa- Similarly, for the second cycle with q = qb, the solution has the form 



fb{t) = abVibit) + I3by2b{t) ■ 



(C2) 



Next we note that the new coefficients ab and (3b are related to those of the previous cycle through 
the relations 



"6 = aa2/ia(vr) +/?ay2a(vr) and 



(C3) 



dt dt 

The new coefficients can thus be considered as a two dimensional vector, and the transformation 
between the coefficients in one cycle and the next cycle is a 2 x 2 matrix. In addition, since we have 
the analytic solution for any given cycle from the previous subsection, the coefficients of this matrix 
are known. Finally, since the equation is symmetric with respect to the midpoint t = 7r/2, and 
since the Wronksian of the original differential equation (j48p is unity, the number of independent 
matrix coefficients is reduced from four to two. The transformation can thus be written in the form 
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where the matrix M (defined in the second equality) depends on the value of qa, and where the 
matrix coefficients are given by 
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cos 
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and 



g = — \/Asin(\/Avr) — qcos^{^/XiT/2) . 
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After N cycles with varying values oi q = qk (where k ranges from 1 to N) , the solution retains the 
general form given above where the coefficients are determined by the product of matrices according 
to 







ao 






Jo. 



N 



where M^^^ = JJ Mfc(g;, 



(C6) 



fc=i 



With this solution, we have transformed the original differential equation (with a random element) 
into a discrete map. Further, the properties of the product matrix M^-'^) determine whether or not 
the solution is unstable and the corresponding growth rate. 
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The evaluation of the growth rate thus requires the repeated multiphcation of 2 x 2 matrices 
with randomly sampled qj values, where the matrix elements are defined through equations (|C4p 
and ()C5p . An estimate for the growth rate can be found by considering each cycle to experience the 
maximum possible growth, given the parameters of that cycle. Under this assumption, the growth 
rate - denoted here as the asymptotic growth rate - is given by 



where Hk = H{qk) is given by equation (jBlOp evaluated at q = qu- It is understood that if \Hk\ < 1 
for a particular cycle, then the growth factor is unity for that cycle. Notice that this expression 
reduces to the form 



In other words, the asymptotic growth rate is the expectation value of the particular growth rates 
for the individual cycles. 

In practice, however, the growth rates of the solutions will not be exactly equal to 700- During 
each cycle, the solution must contain an admixture of both the growing solution and the decaying 
solution over that cycle. Further, the solutions must match up at the cycle boundaries. Counter- 
intuitively, this effect tends to make the true growth rates somewhat larger than the asymptotic 
growth rate defined here (see AB07 for greater detail). Numerical experimentation indicates that 
the asymptotic growth rate (as defined in eq. |C7j ) is generally within 20% of the true value, 
although this result varies with the part of parameter space under consideration. The analysis of 
AB07 shows that both the variations in A and the matching conditions across cycle boundaries 
can produce significant contributions to the growth rate. This correction is important for cases in 
which the variance cr^ of the ratios of the matrix elements becomes extremely large, a limit that is 
not expected to apply to astrophysical orbit problems such as those considered here. Specifically, 
this correction to the growth rate has the form A7 = (T^/(87r) in the limit a <^ 1 and the form 
A7 = CooC in the limit cr ^ 1, where the constant Cqo depends on the shape of the distribution 



Most of the discussion in the main text (§5) considers the forcing function Q{t) to have zero 
width (the delta function limit) and to be exactly vr-periodic. Using the development of a discrete 
map found in Appendix C, we can relax the second assumption (the first assumption is generalized 
in §5.2). As before, we consider one cycle at a time, where each cycle has an amplitude qj, but 
the period is now given by fijir. In the development thus far, we have defined the time variable to 
make the period equal to vr (which can always be done as long as the period is constant). In this 
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(AB07). 



Appendix D: Generalization to Aperiodic Variations 
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generalized treatment, we define the time variable so that the mean period is vr. As a result, the 
average of which we denote with angular brackets, is taken to be unity so that {1/fJ-j) = 1. 

The equation of motion now takes the form 

y = 0, (Dl) 

where we have normalized the forcing frequency to have unit amplitude {Q = Q/qj)- Since Q (and 
Q) are vr-periodic, the jth cycle is defined over the time interval < fijt < vr, or < t < 7r//ij. For 
each cycle, we can thus re-scale both the time variable and the "constants" according to 

t fijt , Xj Xj/fij = Xj , and qj qj/fJ-^ = qj , (D2) 

so the equation of motion reduces to the familiar form (see Theorem 1 of AB07) 

y = 0. (D3) 

For example, if we consider the delta function limit, then Q{t) = 5{[t] — 7r/2), and the solution for a 
given cycle is determined by the matrix M(gj), which is defined by equation ()C4p evaluated using 
Xj and (jj to specify the matrix elements. Notice also that we have allowed the parameter Xj to also 
vary from cycle to cycle within this formalism. Since we need to re-scale the value of Aj every cycle 
to account for the period variation (specified by fij), the net cost of evaluation is the same. The 
physical implication of this result is that departures from exact periodicity (in the time interval, 
not the amplitude) lead to an additional effective variation in the forcing strength qj and the base 
frequency (determined by Aj), but do not otherwise affect the instability. In general, an increased 
variance in the distribution of the parameters (Aj, qj) leads to larger growth rates (AB07). 

For the orbits considered in this paper, the period does not vary from cycle to cycle as much 
as the forcing strength (see Fig. [6|). For example, the standard deviation cr^ of the relative (scaled) 
period is generally small compared to unity, so that the parameter fij typically lies in the range 
given by fij = 1 it o"^, where cr^ <C 1. As a result, the corrections to the frequency parameter 
and forcing strength parameter (from eq. |D2j ) have a typical relative size of only 2a ^ <^ 1. The 
variations in qj from cycle to cycle are much larger and dominate the random component of the 
equation of motion for the instability. 



ay 
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Xj + qjQifijt) 
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Aj + qjQ{t) 



Appendix E: Heuristic Discussion of Hill's Equation 

In this Appendix we outline the essential features of Hill's equation, the reasons for instability, 
and the relation of these issues to the solutions found in this paper. In rough terms, the form of 
Hill's equation (or, equivalently, Mathieu's equation) allows for unstable solutions because it acts 
like a driven oscillator. We first note that if the forcing term were absent (Q(t) = 0), then the 
equation would describe simple harmonic motion with a frequency given by \/A. For nonzero forcing 
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functions, one can think of the forcing term as being on the right hand side of the equation, so 
that one obtains a driven oscillator, which can have growing (unstable) behavior. In this case, the 
driving term is complicated by the fact that it depends on the solution, y{t), which in turn depends 
on both the driving term and the natural oscillation frequency \/A- However, this complication 
does not prevent unstable growth. 

Next we note that since there are two frequencies involved, the periodic driving frequency and 
the natural oscillation frequency, the two effects can interfere either constructively or destructively. 
As a rough illustration of this behavior, one can imagine pushing a child on a swing (e.g., Arnold 
1978). If the pushes are in phase with the oscillatory motion, the swing gains amplitude; if the 
pushes are out of phase, then the swing loses amplitude and the oscillations are damped. Because 
of this tuning, or phase requirement. Hill's equation has regions of both stability and instability, 
depending on the relative phases of the two frequencies involved. As result, the plane of parameters 
shows alternating stripes of stability and instability, as shown by Figures [9] and [101 In the delta 
function limit, we can find analytic solutions for a given set of parameters (Appendix B), or for 
a given cycle in the case of random variations (Appendix C). The corresponding solutions, or the 
matrix elements, can be written in terms of sines and cosines of the angle 6 = \/Avr. Here, \/A is 
the natural oscillation frequency and tt is the period (by construction) of the forcing term. The 
combination of these two quantities thus determines the solutions and the growth rates. As the 
frequency \/A varies (for fixed period vr of the driving term), the system goes in and out phase, 
resulting in the alternating bands of stability and instability shown in the figures. 

The tuning described above depends on the amplitude q of the driving term. Much like for 
the case of pushing a swing, a stronger driving term allows for greater error (leeway) in the phase 
matching necessary to maintain unstable growth. Since the required tuning to obtain instability 
becomes less delicate, more of the available parameter space allows for instability, and the regions 
of stability get smaller with increasing forcing amplitude q (Figs. [9] and [TO]) . Nonetheless, the 
regions of stability never vanish completely. The analysis of Appendices B and C illustrates this 
behavior quantitatively for the simplest case of Hill's equation in the delta function limit. 
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